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ABSTRACT 

We use cosmological adaptive mesh refinement (AMR) code Enzo zoom-in simulations 
to study the long term evolution of the collapsing gas within dark matter halos at 2 ;. 
This direct collapse process is a leading candidate for rapid formation of supermassive 
black hole (SMBH) seeds. To circumvent the Courant condition at small radii, we apply 
the sink particle method, focusing on evolution on scales ^ 0.01 — 10 pc. The collapse 
proceeds in two stages, with the secondary runaway happening within the central 
10 pc. The sink particles form when the collapsing gas requires additional refinement 
of the grid size at the highest refinement level. Their growth is negligible with the 
sole exception of the central seed which grows dramatically to Mgeed - 2 X 10^ Mq 
in ^ 2Myr, confirming the feasibility of this path to the SMBH. The variability of 
angular momentum in the accreted gas results in the formation of two misaligned 
disks. Both disks lie within the Roche limit of the central seed. While the inner disk 
is geometrically thin and weakly asymmetric, the outer disk flares due to turbulent 
motions as a result of the massive inflow along a pair of penetrating filaments. The 
filamentary inflow determines the dominant Fourier modes in this disk — these modes 
have a non-self-gravitational origin. We do not confirm that m = 1 is a dominant 
mode that drives the inflow in the presence of a central massive object. The overall 
configuration appears to be generic, and is expected to form when the central seed 
becomes sufficiently massive. 

Key words: methods: numerical — galaxies: formation — galaxies: high-redshift — 
cosmology: theory — cosmology: dark ages, reionization, first stars 


1 INTRODUCTION 


Supermassive black holes (SMBHs) formed early in the 
history of the universe. Luminous quasar s have be e n ob¬ 
served up to redsh i fts z ^ 7.1 ( e.g., I Fan et aP l2003l : 
iMortlock et al.|[201ll : IWu et al.|[2015l b and have been esti¬ 
mated to host SMBHs in excess of 10® M©, just 700 Myr 
after the Big Bang. Unless the SMBHs are primordial, 
they must have formed by the accumulation of matter 
during the epoch of galaxy format ion — either as rem¬ 
nants of the Population HI stars (e.g Haiman Loebll 200ll: 


lAbel et ^ I 2 OO 2 I : iBromm LarsoiJ12004 ; IVolonteri ReesI 
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I 2 OO 6 I : iLi et al.ll2007l : IPelupessv et ahllioOTl l , or the end prod¬ 
ucts of gas collapse into dark m atter (DM) halos with 
virial tempera tures ^ 10^ K (e.g.. [ Haeh nelt Rees| 19931: 

■ i|2 


Oh Sz HaimanllioO^: Bromm LQebll2003l:IVolonteri Reed 


20061: iBeeelman et al.l I 2 OO 6 : iBeeelman & ShlosmanI 20091: 

Milosavlievic et al.l 20091: Maver et al.l 2010l: Sch 

eicher et al. 

2010 

iHosokawa et al.ll201ll: 1 Johnson et al. I 2 OII 

: IChoi et al. 


for radiative feedback during formation of Pop HI stars low¬ 
ers their main sequence masses to be more in line with those 
of normal stars. The black hole remnants are downsized as 
well, to ~ 10 M©, making it exceedingly difficult to explain 
the growth of early SMBHs from Pop HI seeds. Furthermore, 
recent results from Planck have pushed the beginning of the 
Pop HI epoch forward from ~ 400 Myr to ~ 560 Myr after 
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the Big Bang te.g.- IPIanck Collaboration et a,l.|l2fl1 fili . Taken 
together, these arguments Favor substantially more massive 
SMBH seeds, ~ 10^ — 10^ M©, as the direct collapse models 
suggest. 

Important details of direct collapse remain unclear, e.g., 
how does the process evolve when the flow becomes opti¬ 
cally thick to internally produced radiation? How is angular 
momentum redistributed during the optically-thick regime? 
Does the collapse lead to a hydrostatically supported super- 
mass ive star (SMS) — a precursor to SMBH seed formation 
(e.g., lBegelmai]|201Ql l — or does the flow remain disky and 
by-pass the SMS and the associated thermonuclear stage 
jChoi et, alJbni.'tl ? 


On the other hand, we now understand some of the 
details of the collapse in the optically-thin regime. The frag¬ 
mentation of the flow is damped by the dominant gravi¬ 
tational potential of the host DM halo and by virial su¬ 
personic turbulence — both of which substantially increase 
the Jeans mass of the gas. Furthermore, the DM acts as 
a sink of the angular momentum from the collapsing gas, 
both in isola ted collapse models and in the cosmological 


context (e-g-, Wise et m]| 2008I : iBegelman Shlosmanll^OOl : 


IChoi et al.ll2^.l2015l^ 


Studies of the direct collapse to form the SMBH seeds 
involve computationally-intensive efforts to follow the hy¬ 
drodynamics and the thermal state of the collapsing gas. 
Processes like fragmentation, star formation, redistribution 
of angular momentum in the collapsing flow, development 
of a supersonic turbulence, etc., can be addressed only nu¬ 
merically. 


Typically, the collapse region extends over many orders 
of magnitude in radius, from a few kpc down to ~ 1 — 100 AU. 
Even sophisticated methods face difficulties in dealing with 
such an extended dynamic range. As the collapse progresses, 
the simulation timesteps become increasingly short, making 
it prohibitively expensive to follow the evolution. Moreover, 
the buildup of the optical depth in the flow can fully or 
partially trap the radiation produced in situ, affecting the 
dynamics and thermodynamics of the collapsing gas. 

Two options have been proposed to overcome these ob¬ 
stacles. The first involves (consecutive) zoom-in simulations 
typically associated with calculations of cosm ological struc¬ 
ture formation and galaxy evolution (e.g., iBecerra et al.l 
l2015l b This method allows one to increase the mass and 
spatial resolution of the model, especially on the smallest 
scales, to introduce various physical processes, on different 
scales, and to avoid artefacts associated with a small num¬ 
ber of particles or grid cells at small radii. iBecerra et al.l 
il2nifJi have assumed that the radiative cooling becomes inef¬ 
ficient at some presumed radius. Such an exponential cutoff 
in the cooling rate naturally leads to a radiation pressure- 
dominated entity — an SMS. By itself this does not re¬ 
solve the issue of radiation trapping in the flow — in a 3-D 
rotation-dominated inflow the radiation, in principle, can 
escape. It does not allow one to study the optically thin-to- 
thick transition in the flow and the nature of the flow inside 
the optically-thick region. Merely increasing the mass and 
spatial resolutions will not resolve the problem and will not 
allow one to prolong the simulation without also introducing 
additional physical processes such as radiative transfer and 
associated dynamical effects. At present time, sufficiently 


powerful radiation hydrodynamics codes are not generally 
available. 

The second option involves the sink particle method, 
successfully applied to star formation and other problems 
to focus on a specific range in radius, preventing the smaller 
scales from forcing an increasingly small timestep, and thus 
enabling one to fo llow the long-term evolution of the col- 


lapsing flow (e. 

g.. iBate et al.l Il995l: Krumholz et al. 

2004; 

Eederrath et al. 

boidilWang et al.ll2010l: Tevssier et al.l 

2011 ; 

Gong & Ostrikerll2013lb In this method, gas on small spatial 


scales is replaced by the sink particles, which can grow in 
mass as the surrounding gas evolves. 


In this paper, we apply the sink particle method to the 
direct collapse problem in a cosmological framework. We 
focus on the long-term evolution of the flow and on associ¬ 
ated processes on scales of ~ 0.01 — 10 pc, deep inside the 
DM halo gravitational potential. The choice of this scale is 
based on purely numerical reasons, in order to circumvent 
the exceedingly small timestep. It also allows us to avoid 
additional physical processes on the small scale, such as ra¬ 
diative transfer, without any loss of generality. 

S tarting with our c osmological simulations of direct col¬ 
lapse (|Choi et al.ll2015l l. we create sink particles in the flow 
(as described in Section[2]), and analyze motions and growth 
of these particles and their effect on the collapse itself. Choi 
et al. has followed the collapse down to ~ 10“^ pc scales, 
where one anticipates trapping of the continuum photons 
produced in the flow and the formation of a photosphere. 
Of course, the Lyman a photons are trapped or at least par¬ 
tially trapped already at larger radii. Here we ignore scales 
smaller than ~ 0.01 pc, which allows us to continue follow¬ 
ing the collapse on scales > 0.01 pc, over timescales much 
longer than those simulated by Choi et al. 

This paper is structured as follows. Section[2] provides 
the details of our numerical technique and explains the sink 
particle method implemented here. We present our results 
in Section[3] and discuss them in the last section. 


2 NUMERICAL TECHNIQUES 
2.1 Simulations 

We use the Eulerian adaptive mesh refinement (AMR) code 
Enzo-2.3, whic h has been tested exten s ively and is pub¬ 
licly available (iBrvan NormanI Il997l : iNorman Brvai] 
1 19991 : iBrvan et al.l l2014r ). Enzo uses a particle-mesh N- 
body method to calculate the gravitational dynamics, in¬ 
cluding collisionless DM particl es, and a second-or der piece- 
wise parabolic method (PPM, IBrvan et al.l Il995l i to solve 
hydrodynamics. The structured AMR used in Enzo places 
no fundamental restrictions on the number of rectangular 
grids used to cover some region of space at a given level 
of rehnement, or on t he number of levels of refinement 
([Berger Colellalll989l b A region of the simulation grid is 
refined by a factor of 2 in lengthscale, if either the gas or DM 
density become greater than po,gas,dmA^, where po,gas,dm is 
the cosmic mean density for the gas or DM respectively, 
A = 2 is the refinement factor, and I is the maximal AMR 
refinement level. 

The Jeans length has be en resolved by at leas t 16 cells 
in the simulations. Hence, the lTruelove et all (|l997l i require¬ 
ment for resolution of the Jeans length (see Section 2.3), 
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i.e., at least four cells, has been superseded substantial!; 


(e.g., Sur et alfeoid : iFederrath et aljfioill : iTurk et ^l201 

iLatif et al.ll2013h . 


i 


Enzo follows the non-equilibrium evolution of six 
species: H, He, He+, He++, and e~ (|Abel et al.lll997l : 
lAnninos et alJ Il997l l in a gas with a primordial composi¬ 
tion. It calculates radiative heating and cooling following 
atomic line excitation, recombination, collisional excitation 
and free-free transitions. Radiative losses from atomic cool¬ 
ing are computed in the optically-thin limit. 


Several radiation processes which prevent H 2 forma- 


tion in the collapsing 

: flow have been proposed (e.g.JOmukai 

2001 

: Soaans & Silkll 

20061: Schleicher et al. 

2 OIOI: Latif et al. 

2011 

: Choi et al.ll2013l: llnavoshi et al.ll2014 

: ISugimura et al. 

2014 

). In this work 

we assume that H 2 does not form and 


exclude the chemistry and cooling related to H 2 , which sim¬ 
plifies the chemical evolution. 


2.2 Zoom-in initial conditions 


the gas decouples, its central density increases rapidly. Be¬ 
low we describe the specifics of this method. 


Various works have i nvoked sink particles in order to 


advance their goals fe.g.. iBate et al.l 19951: Krumholz et al. 

2 oo4 

Eederrath et al.ll2010l:IWang et al.l 2010l:lTevssier et al. 

201 1 |; 

Gong & Ostrikerl 2013lj. We implement the sink parti- 


cle method largely based on IWang et al.l (|2010li . who stud¬ 
ied massive star formation using Enzo AMR simulations. 
We require a sink particle to form when the cell violates the 
refinement criterion at the highest refinement level in the 
collapsing part of the flow. 

When this happens, a sink particle is inserted at the 
center of the cell. The initial sink particle mass is computed 
based on the mass exceeding the maximum allowed density 
at a maximum refinement level. In other words, each cell 
density has a maximal value that does not violate the refine¬ 
ment criterion. The initial sink particle velocity is calculated 
based on linear momentum conservation. 

After a sink particle is formed, it accretes the gas from 
its host cell according to following prescription: 


We follow the long-term dynamical evolution of the collaps¬ 
ing gas within a DM halo in a fully cosmological environment 
and subject to atomic cooling. To s atisfy the resolution re¬ 
quirement, we use the MUSIC code (|Hahn V Abell[201l| j to 
generate the cosmological zoom-in initial conditions (ICs). 
MUSIC uses a real-space convolution approach in conjunc¬ 
tion with an adaptive multi-grid Poisson solver to generate 
highly accurate nested density, particle displacement, and 
velocity fields suitable for multi-scale zoom-in simulations 
of structure formation in the universe. 

Generating a set of “zoom-in” ICs is a two-step process. 
Eirst, we generate lh~^ Mpc comoving 128^ DM-only ICs for 
the pathfinder simulation and r un it without AMR until z = 
10. Using the HOP group finder (|Eisenstein V Hul]ll998l i. we 
select an appropriate DM halo, whose mass is ~ 10^h~^ M© 
at 2 : = 10. Secondly, we generate lh~^ Mpc ICs with 512^ 
resolution in DM and gas in the zoom-in region. Since we 
use the same random seeds of these ICs as the ICs at the 
first step, the phases of both ICs are identical. The zoom-in 
region is centered on the selected halo position and is set to 
be large enough to cover the initial positions of all selected 
halo particles. 

The ICs are generated using WMAP5 cosmology: Qa = 
0.721, Qm = 0.279, Ub = 0.0445, h = 0.701, cr = 0.807, 
and Tie = 0.961. In the following we use R for spherical 
coordinates and r for cylindrical ones. 


2.3 Sink particle algorithm 

Sink particles are used in order to weaken the Courant con¬ 
ditions restrictions on the timestep. The restriction can re¬ 
sult from the increase in the local density of the gas and 
demand for additional refinement (i.e. spatial resolution). 
Under these circumstances, it is nearly impossible to study 
the large-scale, long-term evolution of the direct collapse 
within DM halos. To circumvent these numerical obstacles, 
we adopt the sink particle method to model the gravita¬ 
tional collapse. This is especially relevant for the present 
work, which aims at long-term evolution of the direct col¬ 
lapse. Before the gas density exceeds that of the DM density, 
at t ~ 350 Myr, sink particles are not used. However, after 


Msp — 47rpout^B^l-2544c^g^^ ( 1 ) 

where Msp is the mass growth rate of a sink particle, 
Pout = Pceiimin[1.0, (/ceii/r’s)^'^], /cell is the the cell size, 
tb = GMsp/(Cceii + is the Bondi accretion radius, Cceii 
is the sound speed in the cell, and I’rei is the relative velocity 
between the host cell and the sink particle. 

Merging between sink particles is a necessary process in 
order to accurately estimate their mass growth rate and to 
reduce the computational cost. In principle, the sink parti¬ 
cles represent overdense gas clumps in our simulations whose 
internal evolution we ignore. We require two sink particles 
to merge when their separation becomes smaller than 5/ceii, 
and assume that the less massive sink particle merges with 
the more massive one. The merging process is quite insensi¬ 
tive to the definition of the cr itical separation between par¬ 
ticles fe.g.. IWang et al1l2010l j. We have confirmed this by 
comparing various maximal refinement levels. 

In this paper, we use three different maximum refine¬ 
ment levels that affect the overall simulation resolution and 
the scale of the sink particle formation. The maximal AMR 
refinement level I (see Section iTT]) is 10 in run LIO, I = 12 in 
run LI2, and / = 15 in run LI5. The maximal spatial reso¬ 
lution corresponds to ~ 0.1 pc for run L12 and to ~ 0.01 pc 
for run L15. 

All three simulations hav e been restarted from the 
cosmological run described in IChoi et ^ (|2015l i. at t ~ 
355 Myr, when the refinement level in the gas is / = 8, and 
when the gas-to-DM density ratio in the central ~ 10 pc is 
about unity. They are identical simulations except for the 
maximal AMR level. While our analysis here is based on 
run L12, we use runs LIO and L15 to test convergence of the 
observed evolution. 


3 RESULTS 

The DM halo under investigation has been selected from a 
computational box containing several such objects at z = 10. 
This halo has been re-simulated at high resolution down 
to z ~ 12, or t ~ 355 Myr from the Big Bang, when it 
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reached the virial mass and radius of Mh ~ 2 x 10^h~^ Mq 
and i^vir ~ 10/i“^kpc in comoving coordinates, with a 
cosmological spin parameter A ~ 0.03. At this time, the 
DM densi ty profile can be ap proximated by an NFW pro¬ 
file (e.g., iNavarro et al.l [l9^ b with a characteristic scale 
radius of Rs ~ 30 pc, in physical coordinates, and with 
a concentration parameter c = R^ir/Rs ~ 25. At about 
t ~ 355 Myr, the gas-to-DM density ratio exceeds unity in¬ 
terior to ~ i^s- The gas density quickly establishes a density 
profile slightly steeper than ~ R~^, down to the limiting 
radius of R ~ 10“^ pc, imposing a tight Courant condition 
on the timestep. 

As shown by IChoi etlil ] (l2ni.l l201,^ . direct collapse, 
in both the isolated case and in a cosmological framework, 
can be divided into two parts: the first stage ends when the 
central gas density reaches that of the DM, and the second 
stage begins when the central gas density exceeds that of the 
DM, and the gas effectively decouples from the DM. Dur¬ 
ing this latter stage, an density profile is established all 
the way down to ~ 10“^ pc — the region where the opti¬ 
cal depth to continuum radiation produced by the accretion 
flow becomes larger than unity, the so-called thin-to-thick 
transition. As such a large dynamic range imposes a very 
strict condition on the timestep, and because we focus on 
the optically-thin part of the flow, we introduce sink parti¬ 
cles based on the criteria discussed in Section 2.3. By doing 
this, we remove the smallest spatial scales, r < 0.1 pc (for 
run L12) and < 0.01 pc (for run L15), from consideration. 
Hence, these are the limiting resolutions of the models. On 
the other hand, we gain the abi lity to follow the collapse 
well beyond the time achieved in I Choi et ell (|2015l i . 

3.1 Rapid growth of the central sink particle 

The masses of sink particles introduced in the run L12 cor¬ 
respond to the gas mass of the smallest grid cells. Runs 
LIO and L15, in this respect, are very similar, but their 
resolutions are ~ 1 pc and ~ 0.01 pc respectively. So, the 
newly born particles provide only a negligible effect on the 
flow when they are introduced. They can grow by accretion 
and merging and, at least in principle, their influence grows. 
We find that sink particles appear always within the cen¬ 
tral ~ 1 — 2 pc, measured from the densest gas cell. Their 
number never substantially exceeds ~ 10 in the collaps¬ 
ing flow. They grow, but their individual masses stay be¬ 
low ~ 10^ Mq. The sole exception is the sink particle that 
lies at the center. This particle grows fast, and in less than 
~ 2Myrs surpasses Mgeed ~ 2 x 10® M©, as shown in Fig- 
ure[Tl (We term this sink particle as a ‘seed’ and conjecture 
that its evolution is closely related to the future formation 
of the SMBH seed.) After this time, its growth rate satu¬ 
rates substantially. We separate these phases into dynamic 
and secular growth stages of the central seed. 

For comparison, we also show the evolution of the cen¬ 
tral sink particle in runs LIO and L15. In run LIO, the central 
sink particle starts to grow slightly earlier, but very quickly 
its growth curve merges with that of LI2, even before the 
end of the dynamic stage. During the secular stage, the LIO 
growth rate initially exceeds that for L12. Asymptotically, 
the LIO and L12 growth rates are nearly identical, but the 
final seed masses differ by about 20%. 

Run L15 displays a nearly identical growth for the cen- 



Figure 1. Evolution of the central sink particle mass, Mgeed ? 
forming at the position of the density peak for three different 
AMR runs, LIO, L12 and L15. While LIO exhibits a somewhat 
slower or faster growth of the central seed at various times, L12 
and L15 show virtually no difference. Note, that the initial growth 
of the seed mass is very rapid, so that it reaches Mgeed ~ 2 x 
10® M© in ^ 2 Myr. 


tral particle, except for very early times — growth is trig¬ 
gered slightly later but the rate is higher. We therefore, 
conclude that the resolution of run L12 is sufficient for our 
purpose, because the growth rates in L12 and L15 appear 
identical. 

So, irrespective of the resolution we choose for the sink 
particles, the central seed grows in a very short time to ex¬ 
ceed 10® M©. How is this mass converted into the SMBH 
seed is a subject for future research. 


3.2 The central massive seed and disk formation 

To understand the specifics of the growth of the central seed, 
we plot the evolution of the gas masses within the central 
2 pc and 20 pc, and compare them with the mass of the cen¬ 
tral seed, the total mass of other (off-center) sink particles, 
and the DM within these spherical volumes. Note that the 
sink particles grow only by accreting gas and incorporating 
other sink particles, but not through accretion of DM. The 
total baryonic mass, i.e., the gas and the sink particles, is, 
therefore, conserved. 

In Figure[2^, we observe that the total mass of the off- 
center sink particles in the region remains very small. The 
gas accumulation inside the innermost 2 pc starts around 
358 Myr and peaks at ~ 360 Myr, prior to and concurrently 
with the formation and rapid growth of the central seed. 
After the rapid growth stage, the gas mass inside this re¬ 
gion is remarkably constant at ~ 10^ M©, showing that net 
inflow into this region is negligible. Below, we discuss the 
detailed mass balance in the region. Note, that the gas at 
the peak inflow rate into this region drags in the DM, which 
dominates the gas mass during the secular evolution of the 
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Figure 2. Evolution of masses within spherical volumes with a radius of 2pc (left frame) and 20pc (right frame). Shown are the gas 
mass (blue dashed line), DM mass (black solid line), central seed (dotted line with large dots), and the total of all minor sink particles 
(dotted green line). 


t=357.21Myr 



Figure 3. Evolution of the gas density profile for run L12. The 
density profiles, p ~ exhibit flattening in the central region, 

from a ~ 2 to a ~ 0.5, where the gas assembles into a disk 
structure around the central massive seed. 

central sink. This provides additional stability against the 
gas fragmentation in the region. 

The mass evolution within the central 20 pc paints a 
similar picture, with an important difference (Fig.[2]D). The 
gas mass grows monotonically until the formation of the 
central seed, then drops by an order of magnitude during 
the 2Myr rapid growth of the central seed. This explains 
the source of the mass in this seed. After this time, the gas 
mass continues to accumulate at nearly the original rate. 
The DM here experiences the same adiabatic contraction 



logio R [pc] 

Figure 4. Evolution of the mass accretion rate profile for the 
L12 run. Note the sharp decrease in M from the peak of ~ 0.3 — 
1 Mq yr“^ by about two orders of magnitude at smaller radii, as 
a result of the rotational support achieved near the massive sink 
particle. 

as in FigureE^, and dominates the gas. It is clear, however, 
that the gas will shortly reverse this situation — as shown by 
FigureEJ), about lOMyrs after the formation of the central 
seed, the gas mass in the region will surpass that of the DM, 
yet it will remain well below the mass of the central seed. 
We discuss reasons for the gas evolution in the vicinity of 
the central seed and its global consequences below and in 
SectionUl 

The differences between gas evolution inside the central 
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Figure 5. Evolution of the radial (top) and tangential (bottom) 
velocity profiles, and nt, in the gas for run L12. Note the sharp 
decrease in the central Vr and increase in the vt there — a clear 
sign of disk formation. 




inner region 
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inner region 


inner region 


Figure 6. Schematic picture of various stages in the direct col¬ 
lapse with the rapidly growing central seed. Stage A: Filamentary 
inflow into the DM halo on scale of ~ 1 kpc; Stage B: Gas density 
exceeds the DM density and the secondary collapse is triggered. 
In the case of current simulation setup, the ‘inner region’ indicates 
roughly 20 pc region from the sink; Stage C: Central sink forms 
and rapidly grows to ~ 10® ^o; Stage D: High angular momen¬ 
tum gas accumulates in a system of misaligned disks (shown here 
edge-on), forming due to the fluctuating angular momentum of 
the collapsing gas. The small inner disk is shown horizontal, and 
the outer disk is inclined. The disks are connected by a symmetric 
warp. See also Figure[8]and explanations in the text. 


region and on larger scales is related to the gas kinematics. 
This is seen in FigurefS] where the slope of the gas density 
within the central ~ 10pc, p ~ r”", changes from a ~ 2 
to ~ 0.5 — the gas density profile has flattened. Because 
we assume that the gas is optically- thin, the tempera ture 
profile remains the same as shown in IChoi etld] (|2015l ) . 

Another by-product of the change in the prevailing kine¬ 
matics, is the evolution of the gas accretion rate profile. The 
peak accretion rate of ~ 1 Mq yr“^ inside ~ 10 pc is reached 
between t ~ 360 — 362 Myr (Fig.U]). This is reflected in the 
very fast growth of the central seed during these 2 Myr. Im¬ 
mediately after this, the radial profile of the mass accretion 
rate, M(r), drops dramatically from its peak by about two 
orders of magnitude or more with decreasing radius. That 
is, at r ~ 10 pc, the mass inflow rate stays approximately 
the same, ~ 0.3 Mq yr“^, but it declines sharply with radius 
to ~ 3 X 10“^Moyr“^ at ~ 0.3pc. Note, that substantial 
amount of gas, ~ 6 x 10® M©, resides in the DM halo outside 
20 pc at the time of disk formation, and that the accretion 
rate onto the halo remains high as well. 

This behavior in the gas density profile and its accretion 
rate within the central 20 pc means that the gas accumulates 
there at rapid pace. At around t ~ 360 Myr, this is related 
to the formation and a rapid growth of the central seed. 
However, after t ~ 362 Myr, this is the direct consequence 
of the accumulation of angular momentum and formation of 
a gaseous disk in this region. 

Evolution of the mass accretion rate is clearly reflected 
in the abrupt collapse of the gas and formation of the central 
seed, and can be followed in the top frame of Figure[5] The 
radial velocity at t ~ 360.6 Myr increases to ^ 70kms“^. 
This is well above the virial velocity of the DM halo, by al¬ 
most an order of magnitude, and can only result from the 
gravitational decoupling of the gas from the background DM 
potential. Rotational velocity around 0.3 pc from the central 
seed has reached vt ~ 160kms“^ by t ~ 363.6 Myr — ex¬ 
ceeding the DM virial velocity by a factor of ~ 16 (Fig.[5l 
bottom frame). Hence the depth of the potential well due 
to the massive seed has substantially surpassed that in the 
DM. Basically, this corresponds to the following sequence of 
events: onset of gravitational collapse in the gas on scales 
^ 20 pc, gravitational decoupling of the gas from the DM, 
formation of the massive seed, and continuing accretion onto 
this seed. At the same time, the increased rotational support 
for the gas can be clearly observed in a decreased radial in¬ 
flow velocity to Vr < 10kms“^, and a corresponding sharp 
increase in the tangential velocity, I’t, which develops a Ke- 
plerian profile inside the inner r ~ 50 pc, as seen in Figure[5l 

The central seed has reached Mseed ~ 10^ M© at 
t ~ 360.1 Myr. Strong compression of the collapsing gas, 
shocks, and radial filaments can be observed at this time. 
The disk first becomes visible at ~ 361.5 Myr, when Mseed ~ 
1.6 X 10® M©, and can be seen in two out of three projec¬ 
tion planes, xz (face-on) and xy (edge-on), initially on scales 
of the inner 2 — 3 pc. It grows rapidly in size, as a strong 
gaseous bar with two associated open spiral^] become dom¬ 
inant around ~ 362.1 Myr. By ~ 362.7 Myr, the disk radius 
has reached ~ 5 pc. Around ~ 363.7 Myr, one starts to dis¬ 
tinguish between an inner, geometrically thin disk, 2-3 pc in 


^ Such spirals signal efficient angular momentum loss by the gas. 
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Figure 7. Gas column density normal to main planes at the end of the simulation, t ~ 369 Myr, for run L12. The left frame shows the 
inner edge-on disk (position delineated by the white arrows) embedded in the outer disk which is inclined by about 45° in the yz plane. 
The middle frame shows the same figure with the inner disk being face-on in the xz plane. The right frame shows another projection of 
the inner edge-on disk and the inclined (toward the observer) outer disk in the xy plane. The inner warp connection between the disks 
can be clearly seen in this projection. The position of the central seed is shown by a black point. The filaments which feed the outer disk 
are clearly visible. 


radius, and an outer, much thicker disk, which is misaligned 
by about 45° with inner one. The overall radius of this con¬ 
figuration is about 10-12 pc. Formation of the outer inclined 
disk is the result of gas influx with a different direction of 
angular momentum into the region. 

The stages of direct collapse are depicted schematically 
in Figure[6l Stage A of the collapse corresponds to the ini¬ 
tial inflow into the DM halo. Stage B corresponds to the 
gas density exceeding that of the DM and dynamically de¬ 
coupling from the background DM potential. Stage C culmi¬ 
nates with the formation of the central seed and its runaway 
growth, and the last stage is associated with the formation 
of the misaligned disk system due to the fluctuating angular 
momentum of the accreting gas. 

This configuration, of inner and outer disks being 
strongly misaligned, appears stable for the next ~ 6 Myr 
— a time period which corresponds to about 20-25 rota¬ 
tional periods at ~ 15 pc. By the end of the simulation, 
t ~ 369.12 Myr, the inner disk has a radius of 3-4 pc and 
the outer disk about 10-12pc (Fig.fT]). The thickness of the 
inner disk is ^ 1 pc, and that of the outer disk varies from 
1 pc to ~ 3 pc — this disk flares as seen in Figure[8] (bot¬ 
tom frames). The inner disk is connected to the outer one 
with a clearly visible symmetric, integral-shaped warp. 

We do not observe any fragmentation in the inner and 
outer disks. Why is the fragmentation suppressed in the 
disks, where a high-density, low-temperature gas resides? 
The answer to this question can be found in Figure[2^. After 
the formation of the massive central seed, the gas inflow into 
the central 2 pc is minimal, while M is high into the outer 
disk (see also Fig.|4]). The mass within the region, R ^ 20 pc, 
is dominated by the central seed and by the DM that was 
dragged inward by collapsing gas, in adiabatic contraction. 
The ratio of the gravitational acceleration due to the gas to 
that of the central seed, Mgas(< R)/Mseed is < 10“^ for the 


inner disk and < 10“^ for the outer disk, even at the end of 
the simulation. 

The necessary condition for the gas to collapse in the 
presence of a tidal field is that its density has to exceed 
the mean density associated with the tidal field. In our case 
this leads to the Roche limit radius of ~ 10 pc. Moreover, 
Toomre’s parameter Q :::> 1 in both disks, due to the pres¬ 
ence of the central seed, and the cooling floor temperature 
of atomic gas. The penetrating filaments also inject gas and 
stir turbulent motions in the outer disk (Fig. (9]). So, at least 
for ~ 10 Myr after the formation of the seed, and until the 
gas builds up anew in the region (Fig.[2]D), the fragmentation 
in the disks will be heavily suppressed. 

To visualize the morphology of the gas inflow, we pro¬ 
duce density slices in three projection planes, xy, yz and 
xz, using nested boxes, from 4kpc on a side to 10 pc, at 
t ~ 369 Myr (Fig. (8]). The inner and outer disks are visi¬ 
ble in the last two frames, i.e., 40 pc and 10 pc. The inner 
disk is geometrically-thin and seen edge-on in the xy and 
yz planes and face-on in the xz plane. It is surrounded by 
the geometrically-thicker disk, whose mid-plane is inclined 
to that of the inner disk by about 45° in the yz plane. Spiral 
structure is clearly observed in both disks. The accreting gas 
is shocked at the disk surface, as can be seen in the bottom 
frames of Figure[9] (left and middle col u mns) . Similar be¬ 
havior has been observed bv IChoi et al.l (|2013l ) for isolated 
models of direct collapse. 


3.3 Evolution of angular momentum in the 

accretion flow: understanding the dynamical 
consequences of the massive seed 

An alternative way to view the growing disk in the cen¬ 
ter is to study the evolution of the angular momentum in 
the collapsing flow (Fig.fTOl). We show the specific angular 
momentum profile, jgas{r), at four representative times. For 
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Figure 8. Gas density slices in three projections at the end of the simulation, t ~ 369 Myr, for run L12. The frames spanning 40 pc and 
10 pc (on a side) exhibit prominent disk features with the xy and yz planes displaying nearly edge-on views, and the xz plane showing 
a nearly face-on view of the forming disk around the most massive sink particle. 
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Figure 9. Vorticity magnitude slices of the gas evolution with projected velocity arrows, at the end of the simulation, t ~ 369 Myr, for 
run L12. Velocity arrows confirm the dominant rotation in the xz plane. Arrows in the xy and yz planes show a shock forming at the 
disk surfaces in the gas which is accreted from outside the disk plane. 
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comparison, we also plot the circular angular momentum 
profiles, jc{r). Initially, the gas has jgas ^ Jc, but increases 
its specific angular momentum with time, as more gas starts 
to move inward. For t — 365.9 Myr, the gas inside ~ 10 pc 
radius from the seed exhibits a nearly maximal allowable 
rotation, jgas ^ jc — a clear signature that the disk com¬ 
pletely dominates the kinematics in this region. 

We now turn to the main mechanism of angular mo¬ 
mentum loss by the gas. If angular momentum were con¬ 
served during the collapse, the gas within the dominant DM 
halo potential would be able to colla pse only by a f actor 
of ~ 10 before forming a disk (e.g., IShlosmanI (20131 . and 
refs, therein). However, the shapes of growing DM halos are 
inherently triaxial, as seen in v i rtually all numerical si mula- 
tions (e.g., lAllgood et ahlliooil : iBerentzen et ahlbood i. The 
gravitational torques which result from such halos remove 
angular mome ntum from the gas and allow the collapse to 
proceed fe.g.. [Serentzen k, Shlosmlinl l2006l h This happens 
even in regions where the gas has decoupled from the DM, 
because for non-axisymmetric mass distributions, the exter¬ 
nal mass can impose torques on smaller radii. 

We have verified the importance of gravitational torques 
on large and small sp atial scales in co smological simula¬ 
tion of direct collapse (jChoi et al.ll2015l h On larger scales, 
r ~ 10 — 50 pc, we have extended this analysis for an¬ 
other ~ 10 Myr, in the presence of the massive central seed 
fFig. llli left frame). A careful inspection of this region shows 
that the high amplitude of the m — 2 mode has its origin in 
a pair of dominant filaments, which feed the outer disk at a 
high mass accretion rate. So, indeed, the development of the 
m = 2 mode with some non-negligible amplitude is expected 
here. Prior to the formation of the central seed, the gas on 
these scales responds to the torques from the DM halo, and 
its response is that of a self-gravitating fluid. However, af¬ 
ter t ~ 362 Myr, the radius of influence of the seed is about 
r’infl 43 Mseed,6Wi?,io where Mseed,6 is the seed mass in 
units of 10® M©, and Uvir,io is the DM halo virial velocity 
in units of 10kms“^. Hence, almost immediately after its 
formation, the central seed dominates the dynamics of the 
region hosting the disks. After this time, the gas response 
to the external torques, both hydrodynamical and gravita¬ 
tional, is that of a non-self-gravitating fluid. 

The morphological evolution of the outer disk is, there¬ 
fore, dominated by the mass influx along the filaments. 
Without this perturbation, the disk would become largely 
axisymmetric, but as the deprojected image of the face-on 
outer disk shows (Fig. [11] middle frame), it is dominated by 
spiral modes, from m = 2 to higher harmonics, with substan¬ 
tial amplitudes. The inner disk, within a cylindrical annulus 
0.5 ^ r ^ 2 pc, and thickness Az = 0.5 pc, exhibits a lower 
amplitude m = 2 and 1 responses (Fig.[TTl right frame), and 
is affected by the filaments to a smaller extent. Note that 
the misalignment of the two disks itself can in principle pro¬ 
duce gravitational torques and induce m = 2 density modes, 
albeit heavily diluted by the central seed. 

Thus, the self-gravity of both disks is severely diluted 
by the gravity of the massive seed in the center. The outer 
disk is perturbed by the pair of external filaments which 
maintains the m — 2 symmetry there. The inflow rate that 
reaches the inner disk is much smaller, comes mostly from 
the outer disk (with infall along the rotation axis), and ap¬ 
pears more axially symmetrized. 



logio ^ [pc] logio ^ [pc] 

Figure 10. Evolution of the specific angular momentum profile 
in the gas, j{r), for run L12 at 358.3 Myr (top left), 361.33 Myr 
(top right), 365.9 Myr (bottom left), and 369.0Myr (bottom 
right). The angular momentum is measured in a cylindrical vol¬ 
ume with a thickness of 10 pc, with respect to the rotation axis of 
the inner disk, which happens to nearly coincide with the y-axis 
of the computational box. The dotted lines represent the instan¬ 
taneous Keplerian specific angular momentum profiles. Note that 
the innermost jgas reaches its Keplerian value at later times. 


How stable is the configuration of misaligned disks and 
why did it form in the first place? The answer to the lat¬ 
ter question is that the misaligned disks form in response to 
the variable angular momentum axis in the accretion flow. 
Future simulations will show whether this configuration is 
limited to a pair of disks, or whether additional misaligned 
rotating flows can form at larger radii. The reason for its 
apparent stability is found in the two filaments which feed 
the outer disk. By the end of the simulation, t ~ 369.12 Myr, 
these filaments have been substantially degraded, and so one 
can expect that the disks will re-align themselves. But grav¬ 
itational torques between the disks are severely diluted by 
the central seed, which helps to prolong the life of this con¬ 
figuration. Of course, continuing collapse will add additional 
mass to the region, and one cannot predict a priori what will 
be the direction of its angular momentum vector. 

To summarize, the massive seed fundamentally trans¬ 
forms the character of the accretion flow, to a large extent 
independently of the numerical resolution, once the radius of 
influence of this seed is resolved. Most importantly, as long 
as it dominates gravitationally, the massive seed dilutes the 
self-gravity of the gas, damps fragmentation, decreases sub¬ 
stantially the angular momentum transfer away from the 
gas, and, therefore, leads to the formation of a disk(s) in its 
vicinity. The scale chosen to break the self-similar collapse 
in our runs, given by the condition to create sink particles, 
determines the position of the disky flow at ~ 0.01 — 10 pc. 
Increasing or decreasing this scale would move the disk in¬ 
ward or outward, respectively, but the appearance of such a 
flow is generic. 
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Figure 11. Evolution of the nonaxisymmetric modes on progressively smaller spatial scales — in the region dominated by the filaments, 
and in the outer and inner disks. Left: Evolution of the Fourier gas density modes m = 1, 2, and 3, normalized to the m = 0 mode 
amplitude for run L12 in the penetrating filaments region. The amplitudes are measured within the cylindrical annulus defined by 
r ~ 10 — 50 pc and thickness Az = 10 pc. Middle: Face-on outer disk and inclined inner disk at the end of the simulation, t = 369.12 Myr. 
Note the effect of the outer filaments penetrating the disk and driving m = 2 and higher modes in the form of spirals. Right: Evolution 
of the Fourier gas density modes m = 1, 2, and 3, normalized to the m = 0 mode amplitude for run L12 in the face-on inner disk. The 
amplitudes are measured within the cylindrical annulus defined by r ~ 0.5 — 2 pc and thickness Az = 0.5 pc. 


3.4 Evolution of turbulent motions in the 
accretion flow 

It is more difficult to recognize the two-disk system in the 
vorticity magnitude slices displayed in Figure[9l Neverthe¬ 
less, the face-on and edge-on inner disk appears prominent 
in the yz and xz projection planes in the bottom frames. 
For the edge-on projections of this disk one observes accre¬ 
tion flows that are directed along its rotation axis. These 
flows appear to be highly turbulent. Generally, the vortic¬ 
ity increases toward the center, as the vorticity maps show 
on various spatial scales. The flow is turbulent both away 
from the outer disk and in its vicinity, though closer to the 
outer disk mid-plane the turbulence decays. Further analy¬ 
sis of the developing turbulence in the collapsing flow will 
be presented elsewhere. 


4 DISCUSSION 

We have used cosmological zoom-in simulations with the 
AMR Enzo code to study the long-term evolution of direct 
collapse that can lead to the formation of SMBH seeds at 
high redshifts. To circumvent the Courant condition on the 
timestep we have invoked the sink particle method. Sink par¬ 
ticles are introduced when the refinement criterion has been 
violated and the resolution becomes insufficient at the high¬ 
est refinement level (see Section[2]). The innermost accretion 
flow on scales of 1 pc forms a number of sink particles after 
~ 361.5 Myr. They have been permitted to grow by accre¬ 
tion and merging. By comparing various spatial resolutions, 
we find that allowing for a maximum of 12 refinements, cor¬ 
responding to 0.1 pc resolution, is sufficient to obtain the 
convergence in the growth rate of the central sink particle, 
i.e., this growth rate remains unchanged when the resolu¬ 
tion is increased. Under these conditions, we have been able 
to prolong the evolut ion of the c ollaps ing gas by another 
10 Myr, compared to IChoi et al.l (|2015l i. who stopped the 
simulations at t = 360.13 Myr (z ~ 12) after the Big Bang. 
Our main results show that. 


• The masses of these sink particles have remained below 
~ 10^ Mq, with the exception of the central particle which 
grows rapidly to Mseed ~ 2 x 10® M© in less than ~ 2 Myr, 
due to gas accretion and merging with other particles. 

• This growth coincides with the cessation of gravita¬ 
tional collapse on sub-pc scales and the formation of two 
misaligned gas disks — the inner one after ~ 361.5 Myr, 
and the outer one after ~ 363.7 Myr. 

• Formation of misaligned disks is related to the variabil¬ 
ity of the angular momentum direction in the region. 

• As we discuss below (see also Section lffi^ . the appear¬ 
ance of the massive object (hereafter central seed) in the 
center is responsible for the cutoff of collapse. 

• On larger scales, however, the collapse continues unim¬ 
peded. At this point, we make no assumptions about the 
nature of this massive object, but note that its rapid forma¬ 
tion confirms the feasibility of this path to an SMBH. 

We have analyzed the co-evolution of the flow with the 
growth of the central seed, on scales comparable to the size 
of the forming disk. For this purpose, we have computed the 
evolution of Fourier modes on scales ~ 0.5 — 2 pc in the inner 
disk, 2.5 — 10 pc in the outer disk, and 10 — 50 pc outside the 
disk system. The low m = 1 to 3 mode amplitudes remain 
moderate, in the range of 0.2 — 0.4 (relative to m = 0) in the 
inner disk, and the outer disk is dominated by m = 2 spiral 
modes. They do not show any additional growth which can 
be expected if the fluid is self-gravitating. In both disks, we 
estimate that gas self-gravity is not important in compari¬ 
son with the gravitational acceleration due to the massive 
central seed, and show that the fragmentation in the disks is 
damped due to the sharp increase of the Jeans masses there, 
after the formation of the massive seed fSection l3.2p . In other 
words, the disks lie within the Roche limit of the central 
seed. However, the inflow into the outer disk is dominated 
by a pair of filaments which maintain the m = 2 perturba¬ 
tion for a long time, while shear tends to axisymmetrize it. 
The mass inflow there is high, M ~ 0.3 M© yr“^. The inner 
disk is fed at a much smaller rate, M ^ 10“^ M©yr“^, by 
the outer disk. Under these conditions of severely diluted 
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self-gravity, mode-coupling does not operate and the non- 
axisy mmetric (and mostly m = 2) modes quickly saturate 
(e.g., IChristodoulou et alJIlQQ^ l. 

Because the innermost accretion rate is substantially 
below that on larger spatial scales, the gas will accumulate 
in the disks, as discussed in section[321 and shown in Fig- 
ure[2l We expect that this mass growth can lead to some 
local gravitational instabilities in the disks, but not to the 
resumption of a massive accretion onto the central seed. 

A number of important issues are related to this evo¬ 
lution. Probably the most interesting one is how the col¬ 
lapse proceeds on smaller scales — this has been deliber¬ 
ately ignored in the current work, with the introduction of 
the seed particle algorithm. We anticipate that the cessation 
of the accretion process observed here justifies our conclu¬ 
sion that the characteristic mass of the forming SMBH seed 
is M. ~ 10® M©. 

Another issue is related to the fate of the inclined disks 
configuration, which formed as a result of gas flowing into 
the region with a variable orientation of angular momentum. 
Most probably this configuration will be destroyed by the 
same process that created it in the first place. In any case, 
if the SMBH seed forms before the next spike in the accre¬ 
tion rate, the radiation and/or mechanical feedback could 
already modify substantially the conditions in the inflowing 
gas. 

The presence of the massive central object, can be fa¬ 
vorable for the growth of the m = 1 mode in the surrounding 
accretion disk, when the disk mass is not negligible compared 
to the massive object — hence wh en disk self-gravity is not 
severely diluted by the object f e.g.. I Adams et al1ll989l b Fur¬ 
thermore, for stellar disks, which are dominated by the cen¬ 
tral object and hence reside in the Keplerian potential, the 
m — 1 mode, when triggered, can be l ong-lived, due t o the 
absence of precession in this potential (|Tremain l ll995[l . Un- 
der these conditions, the individual stellar orbits can be 
weakly “glued” by their gravitational attraction and can 
maintain this non-axisymmetric mode. This mode has been 
suggested to be responsible for the ma ss accretion rate onto 
the S MBH under certain conditions (|HoDkins QuataertI 

l20ldl . 

However, one expects that the m = 1 mode will not 
be able to grow when strong time-dependent perturbations 
are present, which destroy the orbital correlations in the 
Keplerian gaseous disk, i.e., in a disk with the ratio of the 
central mass-to-disk mass much greater than unity. We do 
not observe any substantial growth of m = 1 mode and any 
mass accretion flux associated with this mode. 

Fast growth of the central object will have dynamical 
consequences for the collapsing flow. Among these, the most 
important appears to be the inability of the flow to get rid 
of its angular momentum. Gravitational torques from the 
DM serve as the main mechanism extracting angular mo¬ 
mentum from the gas, but close to the massive central seed, 
the torques are diluted by the gravitational monopole and 
become inefficient. The gas acquires rotational support and 
settles in a disk-like configuration around the center. Its frag¬ 
mentation is damped by the tidal forces of the central seed 
in the presence of the floor temperature of atomic gas. 

Damping of the nonaxisymmetric instabilities in the 
disk leads to the inability of the gas to transfer its angu¬ 
lar momentum to the outer gas or to the surrounding DM. 


This results in the accumulation of gas in this region and 
the growth of the disk. In fact, the rapid growth of the cen¬ 
tral object is the prime cause of the formation of the disks 
on these scales and of their lack of fragmentation. After the 
initial period of rapid growth, the massive central object 
continues to grow but at a much slower pace, by a factor of 
~ 2 - 3 (Fig.[T]). 

The corollary of this evolution is the formation of a mas¬ 
sive object in the center of the direct collapse region. The 
resolution scale in our model is chosen somewhat arbitrar¬ 
ily, just to allow for the second stage of the collapse to be 
resolved spatially. But our results bring up an important is¬ 
sue: how does the nature of the accretion flow changes with 
the formation and growth of the central massive object? The 
resulting formation of the gaseous disk around this object 
is a direct consequence of dramatically reduced efficiency of 
angular momentum transfer within the radius of influence 
of the massive center, fSection l3.3D . Of course the com¬ 
pactness parameter of such an object, Mseed/Rseed, where 
Rseed is the seed size, will determine the inner edge of the 
forming disk, but not the position of its outer edge, as long 
as it is resolved numerically. Here Rseed is the size of the 
seed or the associated numerical resolution. 

We have also assumed that the molecular gas, H 2 , is 
not present during direct collapse, and, therefore, the floor 
temperature of the atomic gas lies around a several thou¬ 
sand degrees. We justify this by the presence of the Lyman- 
Werner co ntinuum which originates in nearby ste llar popula¬ 
tions fe.g.. lAgarwal KhochfaJl20I5l j or in situ (|Choi et al.l 
l2013h . 

Another important issue is our neglect of the radiation 
feedback in the region of the disk. To address this, it is cru¬ 
cial to know the nature of the massive object, and the source 
of the radiation field. If this object is an SMS fueled by the 
thermonuclear rea ctions in the cor e and accretion energy in 
the envelope fe.g.. iBegelmanllioiOl j. this will specify the ef¬ 
fective temperature, spectral energy distribution and geom¬ 
etry of the radiation field. Alternatively, if the SMS stage is 
bypassed (see Section I), the rate of locally produced energy 
will be very different and the radiation field could become 
anisotropic, with most of it escaping in the preferred direc¬ 
tion, i.e., the rotation axis. We also note that very large 
accretion rates encountered during direct collapse can make 
the radiation feedback much less efficient dynamically. 
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